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Abstract. We have developed a multi-orbital approach to compute the electronic 
structure of a quantum impurity using the non-crossing approximation. The 
calculation starts with a mean-field evaluation of the system's electronic structure using 
a standard quantum chemistry code. Here we use density functional theory (DFT). 
We transformed the one-electron structure into an impurity Hamiltonian by using 
maximally localized Wannier functions (MLWF). Hence, we have developed a method 
to study the Kondo effect in systems based on an initial one-electron calculation. We 
have applied our methodology to a copper phthalocyanine molecule chemisorbed on 
Ag (100), and we have described its spectral function for three different cases where 
the molecule presents a single spin or two spins with ferro- and anti-ferromagnetic 
exchange couplings. We find that the use of broken-symmetry mean-field theories such 
as Kohn-Sham DFT cannot deal with the complexity of the spin of open-shell molecules 
on metal surfaces and extra modeling is needed. 
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1. Introduction 

Since the study of the Kondo features of Ce ad-atoms on Ag (lll)[lj and Co ad-atoms 
on Au (111) [2], the scanning tunnelling microscope (STM) has become a privileged 
tool in the study of surface Kondo physics. The STM is a non-intrusive probe that 
can address adsorbed objects at very low bias with very small currents. Hence, the 
STM basically explores the equilibrium properties of the adsorbed systems. Besides ad- 
atoms, Kondo physics has been revealed in objects of increasing complexity, from single 
adsorbates [H [21 [3] to large organic molecules [H El El [71 [8] , ordered nanostructures |9] 
and ad-atom-ligand structures (lOl [TT] . 

Several recent reports show that organic molecules display a Kondo state due to 
a spin in their extended vr-orbitals even when they are adsorbed on a metal surface 
[3 [HI [121 [13]. This is somewhat of a surprise for two reasons: {i) some of these molecules 
are magnetic in the gas phase because they have a magnetic atom, {ii) an extended vr- 
orbital is expected to have a large overlap with the metal surface likely driving the 
vr-system into a mixed valence regime instead of a Kondo one. Hence, the appearance 
of Kondo physics in these systems depends on a series of parameters where common 
wisdom is likely to fail. It is then of great interest to perform calculations to rationalize 
the particular features of adsorbed large molecules that depend as little as possible on 
adjustable parameters. 

We have implemented an impurity solver for the Anderson Hamiltonian ^14j that 
reads the one-electron structure from a density functional theory (DFT) calculation 
for a given spin configuration of the impurity (in the present case the molecular 
orbitals involved in the spin configuration), uses the DFT hybridization to compute 
the dynamical electron exchange with the substrate, and assumes a single electron 
fluctuation which corresponds to the f/ — oo limit of the Anderson Hamiltonian. The 
impurity solver uses the non-crossing approximation NCA [151 [El [13, HE] in the multi- 
orbital formalism by Kuramoto [16] . The multi-orbital aspects of the Kondo problem 
must be correctly taken into account in realistic accounts of the Kondo problem as 
shown by Kroha and collaborators [191 [20] • 

The NCA is a reliable approximation [21] away from very low temperatures where 
it fails to reproduce the Fermi-liquid behavior in fully screened Kondo systems [22] 
and where spurious spectral feature appears at the Fermi energy [221 123]. At typical 
experimental temperatures NCA is a good choice for this type of calculation because 
it retains all the electronic structure of the one-electron part of the Hamiltonian while 
keeping a correct description of the main Kondo features. In fact, the main limitation of 
our theory rather comes from the use of the customary local or semi-local approximations 
to DFT. Indeed, DFT calculations on electronic gaps lead to discrepancies of a factor 2 
off the experimental gap [21] . This is definitely a big drawback when evaluating Kondo 
temperatures (T^) because they depend exponentially on the molecular level value. The 
Kondo temperature depends on the ratio of the molecular eigenvalue to its broadening, 
due to the molecular hybridization with the metal substrate. We will use this ratio as 
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a parameter. Yet, our procedure retains the symmetry and the relative strengths of the 
one-electron Hamiltonian as given by DFT. A second important approximation of our 
work is that of the infinite intramolecular Coulomb energy, [/ — )■ oo. Extensions of NCA 
to treat finite U [25l|26l|27] lead to considerable improvement of the Tk- However, since 
the LDA based electronic structure impedes the calculation of the T^-, the infinite U 
approximation largely suffices for our purposes. 

In NCA, The Kondo physics is modeled by virtual fluctuations of the impurity 
occupancy. These fluctuations are made possible by coupling the Kondo impurity 
to the substrate. In principle, NCA can treat all possible configurations. The 
f/ — 7- oo approximation is admissible, provided that the self-energy contribution due to 
configurations with two electrons or more is negligible. One can generalize this idea and 
indeed consider two sets of configurations which differ by the addition of one electron, 
energetically separated from other configurations by some large value "t/" [I6]. Hence, 
with NCA we can treat impurities of increasing complexity, where the physics involved 
will correspond to virtual fluctuations among configurations differing in electron number 
by one. 

In this way, Roura Bas and Aligia have treated the singlet-triplet quantum 
transition of an Anderson impurity within NCA [28| 129] . We use a similar approach to 
describe the configurations of a copper phthalocyanine (CuPc) on Ag (100). On this 
surface, CuPc captures one electron from the substrate [30] while maintaining a very 
localized spin on the copper atom at the center of the molecule. Hence, this system is 
properly characterized by a Hamiltonian describing singlet-triplet transitions. In order 
to perform our multi-orbital NCA calculations, we first simulate CuPc/Ag(100) with 
DFT, we transform to a maximally-localized Wannier function (MLWF) basis set [31], 
and solve the DFT Hamiltonian expressed in this basis set with our multi-orbital NCA 
code, selecting molecular orbitals that take part in the Kondo physics. 

The use of MLWF is mandatory to be able to unambiguously transform the DFT 
Hamiltonian into an Anderson-like one as we have shown in [31]. This is perhaps the 
biggest difference with other works using impurity solvers based on DFT calculations 
[32] [33] IM] . Our approach is thus algorithmic, except for the tuning of the molecular- 
orbital levels with respect to their hybridization with the substrate since these are 
quantities where current DFT approaches fail. 

In the following section we give more details regarding the implementation and 
execution of the calculations on CuPc/Ag(100). In section [3l we first show the DFT 
results and their conversion to the MLWF basis set. The MLWF results permit us to 
explore the electronic structure that is involved in Kondo physics. We also evaluate the 
spectral functions for the lowest unoccupied molecular orbital (LUMO) that is actively 
involved in the generation of Kondo spectral features. Finally, we analyze our results 
and conclude this work. 
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2. Method 

The NCA is particularly adequate for describing the electronic structure at different 
energy scales. Our multi-orbital implementation |16] uses the mean-field results of a 
local density approximation (LDA) calculation, transforms the LDA Hamiltonian into 
a MLWF basis which permits us to write an Anderson Hamiltonian, [/ — )■ oo, while 
keeping the full multi-configurational aspects of the problem [31]. The NCA is applied 
on the obtained Anderson Hamiltonian. In this way, the calculated electronic structure 
contains all molecular+substrate information in the presence of Kondo physics. In this 
section, we give details on how this is achieved, with special care in the choice of the 
configurations that will determine the Kondo physics of CuPc on Ag (100). 

2.1. Density functional calculations 

Density functional calculations of gas phase copper phthalocyanine (CuPc) capture the 
relevant electronic and geometric properties of the molecule |35]. Briefly, CuPc is a 
D4/1 molecule, figure [H where the d electrons of the central copper atom are split by 
the ligand field of the surrounding atoms. The ligands capture ~ 2 electrons of the Cu 
atom rendering the molecule in a magnetic d^ configuration [35] . The d manifold is split 
such that the dx2_y2 orbital is singly occupied (SOMO, the singly occupied molecular 
orbital). The following empty orbital (LUMO, lowest unoccupied molecular orbital) has 
TT character and double degeneracy, because it constitutes the Cg representation [35] of 
the point group. 

Previous studies using LDA of CuPc on Ag (100) show that the main electronic 
properties of the adsorbed molecule are retained in the calculation [30]. Here, we take 
on that work and extract the relevant one-electron physics important for determining 
the Kondo state. 

We have used the Siesta code [36], relaxed the molecule and first two surface layers 
to forces below 0.04 eV/A using the geometry of the reference pO]. The calculations are 
periodic, using a super-cell containing five layers of 7 x 7 Ag atoms. This unit cell size is 
converged for computing electronic structure features of the adsorbed CuPc on Ag(lOO) 
|30] . The electronic structure calculations have been done with norm- conserving pseudo- 
potentials [37] and strictly localized, DZP numerical atom-centered basis set, optimized 
for this system [30]. The basis sets are build using the method by Anglada et al. [38] , 
and we have used the improved noble metal surface description using the basis sets of the 
reference [39]. The confinement radii (all in Bohr) of s-type first-(^ basis functions were 
6.06 for hydrogen, 5.20 (C), 5.64 (N), 5.77 (Cu) and 5.55 (Ag). The confinement radius 
of the first-C extended diffuse orbitals was 7.75. Good agreement with our previous work 
[30] is thus attained. 

The Kohn-Sham Bloch functions were calculated on a Monkhorst-Pack grid 2x2x1 
and transformed to the basis of maximally localized Wannier functions [lOj HH W2\ 
(MLWF). Our method for obtaining MLWF from Siesta is described in [31] . 
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2.2. Maximally-localized Wannier functions 

The main reasons for working with a MLWF basis set are: {i) the reduction of the 
computational problem by projecting the M-dimensional Hilbert space of Kohn-Sham 
states, where M is the full dimension of the CuPc plus Ag (100) electronic problem, 
onto a smaller A^- dimensional space which faithfully represents [41] the energy bands 
around the Fermi level taking part in the Kondo effect; [ii) the orthogonality of the 
MLWF basis set that permits us to use standard multi-orbital NCA [HI |22], and {iii) 
the extreme localization of MLWF plus their orthogonality gives us a natural way to 
partition the problem into impurity and substrate subspaces pT] . 

Maximally localized Wannier functions have been used in the study of strongly 
correlated matter (for instance see [l3l HH HH |16]). In some cases, the spread 
minimization is skipped and the Kohn-Sham electronic structure is projected directly 
onto trial orbitals. This strategy has been used in the study of transition metal 
oxides [47]. In our system, we found that the projection onto trial orbitals gives 
unreliable results. The disentanglement method [H] can be contrasted with direct 
selection of certain bands of interest. This is often done in the study of transition metal 
oxides where bands bear strong orbital character [471 HE]- Iii the adsorbate problem, the 
band structure is far more complex, which impedes direct band selection. We conclude 
that the MLWF are an optimal choice for the problem of a large molecule on a metallic 
substrate. 

For the substrate's MLWF, we choose to describe the s — p bands only. In order to 
achieve this, the s —p bands of the slab are generated by interstitial Wannier functions, 
in the same way as the surface state of Cu(lll) was achieved in [31]. This is a reasonable 
selection, as long as Kondo physics is considered, because the d bands start at ~ -3 eV 
below the Fermi energy. This is away from the relevant region in energy with an energy 
scale several orders of magnitude larger than the typical Kondo scale. Additionally, 
(i-states are rather localized, presenting small coupling with molecular states. For these 
two reasons d states can be safely omitted in the description of the substrate electronic 
structure taking place in the Kondo effect. The choice of the interstitial centers for the 
MLWF is delicate because CuPc on Ag (100) displays a C4 symmetry [30] and it is crucial 
to ensure that MLWF do not reduce the symmetry of the problem. An unexpected 
problem of the obtention of MLWF for an extended surface is that convergence to 
MLWF is more difficult as the lateral dimensions of the super-cell are increased. The 
surface had to be repeatedly tested to achieve a good description of its electronic bands 
within 2 eV of the Fermi energy. 

The molecule's MLWF are obtained jointly with the substrate ones by using the 
disentanglement method [4T]. In order to achieve the clear partitioning between molecule 
and substrate, the initial set of trial functions consists of 5 d orbitals of the central copper 
atom (see figure [1]), 16 s states for C-H bonds, 32 pz orbitals for every carbon atom, 24 
orbitals for the nitrogen atoms and 40 s orbitals for C-C and C-N bonds. This Wannier 
function description was tested in the two separate systems: the clean Ag (100) surface 
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Figure 1. Ball-and-stick scheme of a copper phthalocyanine molecule. The central 
atom is the copper one, and the rest of atoms are nitrogen, carbon and hydrogen in 
increasing distance from the central one. The free molecule has a D^h symmetry which 
is reduced to C4 upon adsorption on Ag (100) [30] . 



and the gas phase copper phthalocyanine molecule. 

2.2.1. One-particle Hamiltonian In the reference [31], we showed how to obtain an 
Anderson-like Hamiltonian from a Kohn-Sham Hamiltonian in a MLWF basis set. 
There, the evaluation of the main quantities of an atomic magnetic impurity in a non- 
magnetic host was described. Especial emphasis was put on obtaining the intra-atomic 
Coulomb energy. In the present section, we extend that study to the case of molecules. 
The Hamiltonian terms have to be expressed in a suitable way to be able to apply the 
NCA scheme to account for Kondo physics on molecular orbitals. 

Thanks to the localization and orthogonality of the MLWF basis, the Wannier- 
projected Kohn-Sham Hamiltonian can be organized into four blocks 



where Hmoi contains matrix elements between molecular Wannier functions, the second 
block following the diagonal, Hsiab-, refers to Wannier functions of the silver slab, while 
the off-diagonal blocks, V , involve couplings between molecule and slab. 

The eigenstates of Hmoi are the molecular orbitals now obtained for the MLWF 
transformed Hamiltonian. As we shall see below, CuPc is a magnetic molecule because 
one electron is kept in copper's dx2_y2 orbital. Furthermore, upon adsorption CuPc 
captures one more electron in the first vr-orbital, the doubly degenerated LUMO. Hence, 
this open-shell structure, when hybridized with the substrate, gives rise to the electron 
fluctuations of the Kondo effect. These states are singled out of the problem and we 
define a special subspace for them, J^imp, the "impurity" subspace. The rest of states 
is lumped together into the substrate's subspace, J^subs- 

Hence, we first apply to ([1]) a unitary transformation that diagonalizes Hmoi and 
leaves Hgiab untouched. In this way, we obtain the molecular orbitals of the molecular 




(1) 
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part of the Hamiltonian, so that the first block of ([T]) becomes diagonal. Explicitly, 

TH^^iT^ TV 



'-slab 



(2) 



The unitary matrix T defines the transformation from MLWF of the molecule to 
molecular orbitals. The second step is to choose the molecular orbitals that play a 
role in the Kondo effect. These orbitals are placed in the first rows and columns of ([2]) 
and define the subspace Mlmp- 

the matrix of the Wannier-projected Kohn-Sham 



After this rearrangement 
Hamiltonian reads 



imp 



v 



V 



\ 



H. 



subs 



(3) 



where the first block is diagonal and contains selected eigenenergies of Hmou fhe block 
Hsubs contains Hamiltonian matrix elements in M'subs and the block V contains couplings 
between states in Mlmp and J^subs- The molecular orbitals that do not directly intervene 
in the Kondo effect are included in the substrate. Hence, the full mean-field structure of 
the DFT calculation is preserved within the inner window of the MLWF transformation 
which is 2 eV around the Fermi energy in the present calculation. 

In the present case, the obtention of an Anderson-like Hamiltonian is facilitated by 
the small values of the intrinsic U term in LDA. As shown in [31], one can extract the 
intrinsic U term when writing the full Kohn-Sham Hamiltonian in a Wannier-basis set. 
In that reference, the Co-atom impurity had values of U around 1 eV. For the present 
method to work, one should subtract the intrinsic U term such as is done in the LDA-I-L'^ 
method [49l |50] or in recent parameterizations of model Hamiltonians [51] . Anisimov 
and coworkers show that this intrinsic U term is related to the Hund's rule exchange and, 
thus, is much smaller than the actual Coulomb U values. This is particularly true in 
the case of molecules. In the case of CuPc, we have evaluated the intrinsic [/-term to be 
~ 30 meV, and hence negligible in front of the molecular level values and hybridization 
function. 



2.2.2. Ab-initio calculation of a hybridization function The hybridization function 
has a fundamental role in the impurity physics [52]. It is of uttermost importance to 
do accurate evaluations of this function for numerical applications. In order to achieve 
this, we diagonalize Hgubs, ©, obtaining Bloch states |kn) and Bloch energies ekn, n is 
the band index and k is from the discretized Brillouin zone with 50 x 50 x 1 k-points. 
Convergency on k-points was checked by a four-fold increase of the k-point sampling. 
The couplings of (I3l) are transformed to the |kn) basis accordingly; we label them 
V^„^, where m indexes states in Jifimp- After these manipulations, we can calculate 

Tmm' (W) = ^kn,mKn,m'^i^ " Ckn) (4) 

kn 
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from the one-particle Hamiltonian It is a matrix in J^imp- 

We emphasize that the Hamiltonian (|3]) is obtained by projecting the Kohn-Sham 
electronic structure onto a set of MLWF. The smaller A^-dimensional subspace 
spanned by MLWF is designed in order to represent the selected energy bands around the 
Fermi level, with the same level of accuracy as the original bands of the M-dimensional 
space of Kohn-Sham states. The total number of Wannier functions N is the sum of the 
MLWF of the molecule and of the substrate. The Wannier description has been tested 
in separate CuPc and Ag(lOO) systems. M usually means the highest band output from 
the ab-initio calculation. Interestingly, we found that it is very important to verify the 
convergence in M. Inclusion of bands with energy of even tens of eV above the Fermi 
level is vital and preconditions the correct projection onto the Wannier space. As a 
consequence, the important value of Tmrni^) at the Fermi level turns incorrect if M is 
too small. In the present case, convergence was attained for A^ = 705 and M = 3000. 

Formally, ([3]) becomes the one-particle part of the Anderson Hamiltonian on 
bringing Hgubs to a diagonal form. Although the procedure we have just introduced 
is straightforward and essentially algorithmic, it is not correct in principle, because the 
on-site energies of molecular orbitals in Himp contain certain part of the Coulomb energy 
that is included in the mean-field-like LDA framework. Furthermore, the inadequacy 
of Kohn-Sham energies of orbitals lying close to the Fermi level is well known. In our 
approach, Himp is replaced by a parameterized model Hamiltonian with many-body 
interactions. This is a natural step, in view of the fact that we would not have to deal 
with Hamiltonian partitions like in ([3]) if standard ab-initio calculations included the 
relevant many-body interactions correctly. 

2.3. Impurity electronic configurations 

In this work, the one-particle Hamiltonian Himp is replaced by a Hamiltonian h with 
many-body interactions. This subsection discusses the physical grounds on which h is 
designed. 

The experimental analysis on the Kondo features in CuPc/Ag(100) [1^ indicates 
that the two-fold degenerate Cg LUMO has the main role in Kondo fluctuations. 
However, DFT calculations [30] show that the spin in the SOMO orbital is not quenched 
by charge transfer from the molecule. Hence, the observed Kondo resonance is 
interpreted [I3] as due to a spin S" = 1 of the molecule, formed by an electron in 
LUMO and one electron in SOMO. The LUMO electron is captured from the substrate 
upon adsorption. 

We adopt this interpretation and model the molecule as an impurity with a two- 
electron S = 1 ground state. Furthermore, the lifetime of the SOMO will be assumed an 
order of magnitude longer than the LUMO's, as will be confirmed by ab-initio calculation 
in section 13.1.21 

The electron in the SOMO will be represented by a local magnetic moment which 
interacts via exchange with the hybridized LUMO. The multi-orbital structure of the 
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Figure 2. Considered electronic configurations for the two-fold LUMO {eg): (a) The 
empty orbital (denoted by 6°), (b) singly occupied orbital e^, (c) doubly occupied 
configuration e^, (d) doubly occupied configuration e]^. 



problem appears in the twofold degeneracy of the LUMO. The configurations that we 
will study are then the ones coming from the filling of the LUMO, figure [21 

We choose the relevant configurations by estimating their energetic accessibility in 
the charge fluctuation process. Then, the relevant parameters determining the electronic 
configuration are: the LUMO on-site energy and the charging energy U . The lowest- 
energy molecular configuration for the adsorbed molecule is singly occupied, as suggested 
by DFT calculations [30]. We can estimate the free-molecule U by evaluating the 
energies of the neutral (-Eo), singly {Ei) and doubly negatively charged molecule. 

Assuming a simple impurity Hamiltonian, then the affinity is given by 

Ei~Eo = eL (5) 

and the second affinity by 

Eii -Ej = 2eL + U (6) 

from these equations and the total energy calculations for the free molecular species we 
obtain that U = 2.08 eV. However, U is substantially screened on the metallic surface. 
A constrained DFT calculation for yr-systems on silver surfaces leads to a reduction of 
U to values between 0.5 and 1.0 eV [53]. In the present case we take ei = —0.35 eV 
(see section [3^ . which leads to je^l < U and the doubly occupied configurations e^, 
Cg^, figure |2] (c,d), are higher in energy than the neutral one, figure |2] (a). Hence, we 
will consider — e° fluctuations. 

Hence, the impurity Hamiltonian contains the eigenvalues of the two LUMO 
as given by the diagonalization of the MLWF Hamiltonian ([2]), plus the exchange 
interaction I with the spin of the SOMO: 

h = ei \aa){aa\ — JSi ■ S2. (7) 

a=l,2 a 

We introduce Hubbard operators which automatically restrict the LUMO occupancy to 
zero or one. The sums are over spin and orbital degrees of freedom, the latter indexed 
by a. The second term in the Hamiltonian is the direct exchange interaction involving 
the spin S2 of the SOMO and the spin operator of the LUMO expressed through Pauli 
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matrices r as 

Si = 5^ (^) 5Z l^^'Xa^l • 

era' a=l,2 

If / = the Schrieffer- Wolff projection of the model onto the manifold leads to 
a SU(4) Coqblin-Schrieffer model [5^. The evolution of spectral features with / leads 
to a quantum phase transition of the singlet-triplet impurity problem, that has recently 
received much attention experimentally [SS] and theoretically [2S1 EH]- Here, we will 
study three cases: / < 0, / = and / > 0. 



2.4. Multi-orbital non-crossing approximation 

The above impurity Hamiltonian ([7]) commutes both with the LUMO occupancy 
operator, = J2a kcr)(acr|, and the total spin operator S = Si + S2 showing that it 
properly describes the molecular properties. Due to the hybridization with the substrate, 
the molecule will exchange electrons. Here, we assume just fluctuations between the two 
above configurations, figure [2] (a) and (b) plus the fluctuations of the SOMO spin, as 
included in the definition of h. 

The full Anderson-like Hamiltonian reads 

H = h + J2^^AaC^a + V (8a) 

kcr 

^ = E E (^k.,a4 |0)(acr| + 1^4,, |aa)(0| c^^) . (8b) 

kcr a 

The hybridization part V and the impurity Hamiltonian h are written with Hubbard 
operators which change the state of LUMO from empty to occupied or vice versa. 
Substrate electrons are described by fermionic operators c^^ and energies e^^. For the 
sake of brevity, all substrate-electronic degrees of freedom are encapsulated in the k 
symbol. The band index will be introduced when necessary. 

Let Hq = H — V . Following Bickers |18], we write down a resolvent operator of the 
impurity 

as a result of averaging over eigenstates \ Q) of the non-interacting substrate with energies 
Eq. The substrate partition function is denoted by 2'^"^^. Since R{z) does not mix 
different occupancies, it can be written as a block matrix 



Riz) 



Ri{z) 
R2{z) 



where Ri{z) refers to one-electron occupancy (ie LUMO empty) and R2{.z) acts on two 
electron occupancies of the impurity. 
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These quantities are calculated via self-energies defined by 

Ri{z) ^ 



Ri^z) = . i . . 

Z - h\N=2 - ^2{z) 

Following Kuramoto |TB], we introduce basis states labeled by a for the configurations 
with one electron and /3 for the two-electron configurations. In the Non- Crossing 
Approximation the self-energies for the fixed occupations are given by all diagrams 
without crossings of substrate electron lines, 

13(5' ko- 

Xi?2/3'/3(c^ + ek^)Vka(/3|«) 
a' a kcr 

The hybridization vertex V\s_cr{o:\P) comes when emitting an electron from LUMO to the 
band state ka which is accompanied by impurity transition from /3 to the state with 
label a. Similarly, \4;o-(/3|a) is brought about when annihilating a one-electron state 
a of the impurity and creating a two electron state by absorbing a substrate electron. 
Explicitly, 

VUa\/3) = {a\c^^V\^) (9a) 
VkM(y) = if^lVclja) . (9b) 

It is convenient to re-express the fixed-occupation self-energies (also called bosonic and 
fermionic self-energies in slave-boson approaches ESI |29] ) in terms of a hybridization 
function that is directly related to the level broadening of the impurity configurations 

T{a'(3\f3'a; oo) = Y^ V^.{a'\(3')V^^{/3\a)6{uj - e^J. (10) 

ko- 

Please, notice that we have not included vr or 2tt factors as sometimes is done in the 
literature. In order to evaluate lifetimes a factor 2tt will be added to the diagonal terms 
of the hybridization function (fTOj) . 

We can now represent the fixed-occupation self-energies in the form 

^la'aioo) = Y, I f{oo')T{a'(5\(5'a- J)R2p'p{uj + (11a) 
S2/3/3'(c^) = Y.j f{-oo')T{a'l3\(3'a-uj')Rio.c.'{uj-uj')duj'. (lib) 



This permits us to efficiently perform all computations using fast Fourier transforms. 
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Let us write {3 = {a,S,S^), where a = 1,2 indexes the two orbitals of LUMO, S 
the total spin and 5*^ one of its components. Similarly, a will denote the projection 
of the spin of SOMO on the z-axis, the only degree of freedom of the one-electron 
configurations. 

In the present case, the substrate is non-magnetic, e,^^ and Vko-,a do not depend 
on the electron spin, a, and the expression for the hybridization factorizes into spin 
and orbital parts (see Appendix). This considerably simplifies the equations (|TT|) . 
Introducing 

Ri{uj) = [cu - ^i{cu)]~^ (12a) 

[Rt''{co)]2 = {uj-ea + ^I)Saa' " ^2,aa'{uj) (12b) 
[Rt\^)]al' = i^-^a- \l)Saa' ' ^2,aa'{uj) (12c) 

and defining the orbital part of the hybridization function by (see (lA.ip ) 

r.„'(^) = $^Ka'^ka5(^-ek), (13) 
k 

we can write 



^2,aa'{u)) = j f{-Uj')Taa'{<^')Rl{<^ - w')dw'. 



(14a) 
(14b) 



The resolvent for the two-electron configurations -Rf now depends on the total spin 
and is a matrix in the orbital space. The equations ( fT2lfT4l) constitute a self-consistent 
system. 

The main quantity of interest, the Green's function of LUMO, will be calculated 
from the general time-ordered correlation function of Hubbard operators, defined 
through the thermal average 



G{a'(5\ T\/3a, 0) = -T,( e"^ \a'){/3'\ e""^ 



(15) 



Starting from the latter expression, we average over the singly-occupied configurations 
aa' and take the Fourier transform [16] to obtain the Green's function of LUMO for the 
given spin multiplet 



-Al,ie)R,ie-u)]e-^^de. 



(16) 
(17) 
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The quantities A2^aa' and Ai are the spectral functions of the resolvents and Zi is the 
impurity partition function. We follow Kuramoto [16] and use defect propagators in the 
numerical implementation of ( ITTl) . see Appendix. 

The spectral function of LUMO is calculated by tracing over orbital and spin degrees 
of freedom, 

m = -l^rn J2 [Gfr (c) + 3 G^lfiu)] . (18) 

a 

We see that the general Green's function f|T5|) factorizes into channels of the total 
spin of an incoming electron (hole) and the molecule f fT7|) . This fact has been used to 
model magnetic inelastic effects induced by the STM [561 iS3 158] . 

3. Results: Kondo Physics of CuPc on Ag (100) 

Full detail on the adsorption of CuPc on Ag (100) is given in other publications [301 159]. 
Here, we build on those results and apply our methodology to obtain the electronic 
structure in presence of the Kondo effect. First, we build the Anderson-like Hamiltonian 
following the recipes of section l2.2.1l and next, we solve the multi-orbital NCA equations, 
section 12.41 presenting the spectral function or PDOS onto CuPc doubly degenerated 
LUMO. Special care is given to the evaluation and presentation of the hybridization 
function, Tmm'i^), of section 12.2.21 because most of the relevant Kondo physics is 
associated with the symmetry and values of this function. 

3.1. Hamiltonian for CuPc on Ag (100) 

3.1.1. Molecular orbitals Thanks to the partitioning made possible by the MLWF 
basis, we can extract the molecular Hamiltonian Hmoi ® computed from the LDA 
calculation of CuPc on Ag (100) |30l |59]. Table [T] presents eigenenergies of Hmoi the 
molecular orbitals close to the substrate's Fermi level. The orbitals are then closely 
related to the gas-phase CuPc ones and hence, we denote them by SOMO, LUMO and 
highest occupied molecular orbital (HOMO) labels. For comparison, the results for 
the LDA Siesta calculation of the gas-phase molecule of the same molecular geometry 
are also given. The MLWF qualitatively reproduce the gas phase molecular levels. 
We emphasize that the e (MLWF) refer to molecular orbitals screened by the metal 
but without direct hybridization and the e (Siesta) are for a true gas-phase molecule 
without any screening or coupling with an external metal. 

Figure |3] shows the PDOS onto these four molecular orbitals. The calculation is 
performed for the LDA k-point grid (2x2x1) which is clearly insufficient to account for 
the continuum character of the substrate electronic states. Hence, the peaks have been 
slightly broadened using a Gaussian broadening of 50 meV. Despite the qualitative 
character of this figure, much information can be gleaned from it. The difference in 
widths between SOMO on one hand and LUMO and HOMO on the other hand is 
substantial. The SOMO peak basically presents the numerical width, and is a featureless 
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Figure 3. Density of states projected onto molecular orbitals (PDOS). The projection 
onto the SOMO is given by bold lines, onto the HOMO by dashed lines and onto the 
LUMO by dash-dotted lines. Energies are with respect to the Fermi energy. Majority 
spin is shown in the upper panel and minority in the lower one. The PDOS are 
convoluted with a 50-meV Gaussian. 



peak, while the LUMO shows oscillations and spreads over several hundreds of meV. 
Then, the latter orbitals hybridize more strongly with the substrate. From this figure 
we can conclude that while the LUMO width is in the range of a few hundreds of meV, 
the SOMO is significantly less. This is in agreement with the move involved calculations 
of the hybridization function that we present in the next section. 

The PDOS are in excellent agreement with the PDOS on molecular orbitals from 
DFT results [59]. Both the eigen-energies and the PDOS give support to our method of 
transforming the DFT Hamiltonian to a MLWF and selecting the impurity Hamiltonian 
using (IH]). 

Table 1. Eigenenergies of the molecular orbitals for the impurity Hamiltonian 
evaluated with MLWF using ([5]), compared with the gas phase LDA calculation using 
Siesta, for the same geometry of the CuPc molecule. All energies in the third column 
were shifted, so that the energies of SOMO {\) coincide. Energies of the second 
LUMO state are in parenthesis. The comparison is only indicative that the MLWF 
capture the molecular properties since the two calculations are performed for different 
systems: the MLWF refers to molecular orbitals screened by the metal but without 
direct hybridization and the Siesta column is for a true gas-phase molecule without 
any screening or coupling with an external metal. 



e (MLWF) [eV] e (Siesta) [eV; 



HOMO (t) 


-1.110 


-0.888 


HOMO ii) 


-1.112 


-0.890 


SOMO (t) 


-0.549 


-0.714 


SOMO {{.) 


0.040 


0.040 


LUMO (t) 


0.185 (0.190) 


0.496 (0.498) 


LUMO (;) 


0.211 (0.216) 


0.520 (0.522) 



It is difficult to conclude on the actual occupancies from the PDOS. Indeed, the 



The Kondo signature of CuPc on Ag(lOO) 



15 



PDOS numerical broadening thwarts any precise calculation of orbital occupancies, and 
the definition itself of occupancy of a molecular orbital in a chemisorbed system is 
somewhat arbitrary. As we show in the next section, the actual peak width of SOMO 
is negligible; the SOMO then remains singly occupied as in the gas phase. The LUMO 
peaks in figure [3] cross the Fermi level, hence these orbitals capture charge from the 
substrate. The LUMO are then the only orbitals that participate in charge transfer 
from the surface. 

3.1.2. Hybridization function The lack of continuum in the PDOS calculation is cured 
in the calculation of the hybridization function (jlj. Thanks to the smaller size of the 
MLWF Hamiltonian (|3]), we can now find the Bloch functions for a very dense k-point 
grid. 

The off-diagonal parts of Vmm'{'jj) are very small (< 1 meV). This is a very strong 
result that shows that the substrate does not mix the different molecular orbitals among 
themselves. Hence, every molecular orbital defines an electronic channel of the system. 

The diagonal elements for SOMO and LUMO orbitals are presented in figure HI 
along with the density of states of the substrate. The three curves are very similar in 
shape. By comparing the hybridization functions with the substrate DOS, we notice 
that the hybridization function grows more slowly than the DOS as the electron energy 
increases. This is an effect due to the couplings, V , between the molecular orbitals and 
the substrate electronic structure. However, at higher energy, it is the DOS that controls 
the behavior of the hybridization function with energy. Figure H] shows that the SOMO 
width has to be multiplied by 20 to be comparable to the LUMO one. The SOMO 
orbital has then a very small mixing with the substrate as compared to the LUMO. 
Finally, in the [/ = picture, the FWHM for the LUMO is 2'KTaa{i^ = cl) = 440 
meV and is comparable to the FWHM of the PDOS peak, figure [3l obtained with an 
insufficient k-point sampling. 

We have also evaluated the hybridization function of the LUMO with d orbitals 
by preparing a Wannier basis-set with d orbitals. In the region of interest here (some 
2 eV about the Fermi energy), this hybridization is strictly zero due to the lack of d 
states at these energies. However, when resonant with Ag d-band (3 eV below the Fermi 
energy), the hybridization function becomes larger than the corresponding values for the 
sp Wannier functions. At 4 eV below the Fermi energy, the (i-electron hybridization has 
maximum of 0.032 eV, while the sp one is 0.035 eV. The consequence of this is that 
the spectral function at -4 eV will not be just a simple Lorentzian tail. However, the 
LUMO orbital is some hundreds of meV away from the Fermi energy, and the effect of 
the (i-band contribution to the overall shape of the LUMO spectral function is negligible 
both for Kondo physics and for the one-electron spectral shape. 

Finally, we comment on the problem of broken spin symmetry in DFT. Spin- 
polarized LDA implies two subsystems: minority and majority spin. This in turn 
says that we have two sets of Wannier functions, two distinct substrates, hybridization 
functions, etc. The main effect of breaking the spin invariance in our DFT calculation 
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Figure 4. Diagonal elements of the hybridization function matrix (in eV) as a function 
of the electron energy, w, with respect to the Fermi level. Only the LUMO and the 
SOMO (note the factor 20) are considered. For comparison, the substrate-projected 
density of states is given in arbitrary units. The delta function of ([4]) has been replaced 
by a 5-meV wide Gaussian; the k-point sampling is 50 x 50 x 1. 



is that the SOMO occupancy is nsoMO ~ 1- As a secondary effect, the substrate and 
molecule become slightly spin-polarized as well. However, this effect is perturbational 
[H] and is not an intrinsic property of the bare substrate and impurity of the Anderson 
model. In what follows we drop the minority spin data of the hybridization function 
and restore the spin symmetry. 

3.2. Multi-orhital NCA results 

Our LDA calculations yield a hybridization function for the SOMO, F^, ten times smaller 
than the one for the LUMO levels. From these values we obtain F^ 4.5 meV. We use 
the standard expression for the Kondo temperature [60] 



and take the ab-initio value for the on-site energy from the tabled! The bandwidth is 2D, 
where a rectangular DOS is typically assumed. Here we have taken D as 10 eV because 
the DOS of Fig. 3 integrates to 703 states, while the calculated DOS at the Fermi energy 
is 32.28 eV-^ Hence 2D=703/32.28=21.8 eV, and D turns out roughly 10 eV. We get 
Tk^s of the range 10~^^ eV. Thus, the energy scale that would correspond to a Kondo 
effect on the SOMO orbital (without exchange coupling to LUMO) is unobservable. 

These arguments show that the Kondo coupling in this system is indeed given by 
virtual charge fluctuations of the two-fold degenerate LUMO. Hence, we are dealing with 
the impurity problem described in section 12. 3[ for which we can calculate the spectral 
function according to the section 12.41 

However, the on-site energies of LUMO in (jT]) as given by their LDA values in the 
table [1] lie very close to the Fermi level, which would correspond to a fluctuating valence 




(19) 
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Figure 5. LUMO spectral function for tiie / = case. Tlie wide peak corresponds to 
a charge excitation and the narrow peak is the Kondo peak (shown also in the inset). 

regime. That has not been observed in the experiment [13]. This is related to the fact 
that the LDA energies of LUMO do not reflect the considerable Coulomb repulsion in 
their spin splitting. We fix these deficiencies by rescaling the eL/Vaaif^j) ratio in order to 
achieve a SU(4) Kondo temperature, Tk,l = Dexp {—\eL\/4raa), of ~20 K as observed 
in the measurements [13]. The results presented in this section are calculated with 
= —0.35 eV and Taai^L) rescaled to 0.01 eV. 

Below Tk,l, the model exhibits a rich variety of physics as the value of / changes. 
For / > but smaller than the Kondo temperature Tk^l, the SOMO is effectively 
decoupled from the LUMO and we recover the case of two degenerated LUMO on the 
same footing as the electron spin: we have an SU(4) system as described above. For 
negative /, |/| < the SOMO becomes screened by a two-stage Kondo effect at 

very low temperatures. For positive and large I S> Taa,TK,L, the Kondo physics is 
that of the under-screened Kondo effect, because of the S = 1 molecular ground state. 
Hence, the model will show singular Fermi liquid characteristics [611 [62l [63] in the low- 
energy domain. Finally, in the limit / — > — oo, we obtain a spin zero molecule with 
orbital pseudo-spin, which is over-screened by two substrate channels. The cross-over 
temperature is, however, exponentially smaller than Tk,l, due to the reduced degeneracy. 
Here we present NCA spectral functions in the / = case and in the intermediate 
regimes Taa > l-^l > Tk,l which are dominated by inelastic spin transitions. 

Figure |5] shows the spectral function of LUMO at T = 7K when the exchange 
interaction between SOMO and LUMO is turned off, / = 0. The ground state of h 
is fourfold degenerate in this case. Since each of the two orbitals of LUMO defines an 
independent scattering channel in the substrate (ie Taa'i^jj) is diagonal), the LUMO is 
subject to a SU{A) Kondo effect. The spectral function shows a broad charge-excitation 
peak of the FWHM given by 4 ■ 27rrj^j;^(ei) and a narrow Kondo peak. 

When the SOMO and LUMO spins are subject to a ferromagnetic interaction 
(/ > 0) the ground state of the molecule without couplings to the substrate is an 
orbitally degenerate S = 1. The excited state is a = orbital doublet. Now, two 
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Figure 6. LUMO spectral function for / = 25 meV (bold line). The spin zero and 
spin one channel contributions are presented. The positive-energy satellite corresponds 
to injection of one electron in the excited electronic structure of the system: the spin 
zero channel. The negative-energy peak corresponds to removing one electron from 
the molecular ground state in the S — 1 channel. 

satellites develop at ±/ from the Kondo peak, figure O These satellites are inelastic 
replicas of the Kondo peak since the spin excitation energy is exactly /. Decomposition 
into spin channels yields that the Stokes (+/) satellite is in the 5 = channel, while the 
anti-Stokes peak as well as the Kondo peak are in the S = 1 channel. This result can 
be understood by recalling the meaning of the spectral function. The spectral function 
at T = yields the probability density to inject one electron in the system when a; > 
or to inject a hole when a; < 0. Hence, the positive energy satellite corresponds to an 
excited Kondo effect triggered by the injection of an electron to the S" = state, which is 
an excited state of h. Similar results have been found by Roura Bas and Aligia [281 129] • 

When the interaction is anti- ferromagnetic (/ < 0), the ground state of h is an 
orbital doublet and is followed by six excited states of spin one. A noteworthy feature 
of the spectral function are the steps typical for an inelastic spin-flip transition (see for 
example [561 EZ] ^"^^ references therein). The inelastic steps enhanced by Kondo effect 
have already been explored in the case of singlet-triplet transitions in nanotubes [64] . 
The spin channel analysis yields that the Stokes peak is now S = 1 which again can be 
rationalized by noticing that it corresponds to injecting an electron in an excited S = 1 
state and the transition is enhanced by an excited Kondo effect. By the same token, 
the anti-Stokes peak is S" = since a hole is created in this channel. 

The spectral function, figure [71 shows a small peak on the Fermi level. This peak 
cannot correspond to a spin-flip Kondo effect. The orbital Kondo resonance is ruled 
out by its exponentially suppressed energy scale oc exp (— |eL|/2raa) as compared to 
Tk,l- NCA is known to produce spurious peaks at the Fermi level [151 ESI HH]- The 
temperature scale at which they appear is given by [60] Tg = TK,L/'^[rK,L/'^aaf^^ = 0.18 
K (strictly valid for / = only). All calculations shown here are performed at 7 K, well 
above the pathology temperature T5. We conjecture that the zero energy structure is 
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Figure 7. LUMO spectral function for I = —25 meV. The spin zero and spin one 
components are presented. The two excitation steps correspond to inelastic transitions 
enhanced by the Kondo effect. As in the S = \ case, / = 25 meV, the steps belong to 
a channel with well defined spin. The positive-energy step is given by the spin triplet 
(5 = 1) channel and the low energy peak is due to the S — Q channel. 



related to the problems of NCA to reproduce spin-split Kondo peaks, as detailed in [22] , 

p. 25. 

A common feature of both J > and J < spectral functions (figures 161171) is a 
certain asymmetry of the Stokes and anti-Stokes peaks. This is caused by asymmetry 
of the hybridization function with respect to the Fermi energy, or by the U = oo 
approximation. We have explicitly verified that the hybridization function has very 
little effect on the respective heights of both satellites. Indeed, the asymmetry comes 
from the fact that system is out of the particle-hole symmetry. Then we assign the 
asymmetry of the satellites to the infinite Coulomb repulsion. 

4. Discussion 

The evaluated LDA spectral function, figure IH seems to suggest a fluctuating valence 
state while a Kondo peak have been experimentally found [I3]. Moreover, if we added 
a strong Coulomb term to the LUMO, the corresponding Anderson Hamiltonian would 
correspond to an empty-impurity regime according to the scaling theory [65]. For these 
reasons, we had to rescale the Kondo temperature in our model by shifting the LUMO 
level and changing the hybridization strength as we showed in the previous section. 

These results show that LDA is unreliable to furnish quantitative data. Any ah- 
initio method is bound to failure when trying to estimate the Kondo temperature given 
the exponential dependence on the main ab-initio ingredients: the level position and 
hybridization. Nevertheless, ab-initio calculations can give valuable qualitative input 
since the correct electronic symmetry addressing both the spin and orbital channels can 
be directly obtained from DFT. Our analysis of the hybridization function has yielded 
important information: (z) the existence of two well-defined orbital channels originating 
in the hybridization of the two LUMO with the substrate has been proved because 
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non-diagonal terms of Taa'{^) are vanishingly small, {ii) both LUMO and SOMO are 
partially charged in the adsorbed system, the LUMO due the sizeable values of Vaa that 
leads to charge transfer from the substrate, and the SOMO for the vanishing value of 
T s that leaves its spin unperturbed. Hence, LDA reveals the orbital SU{2) symmetry 
associated with LUMO electron and its substrate channels. 

The failure of the present LDA calculations to yield a realistic ei/^aa ratio can 
be traced back to the complete failure in giving the qualitative Kondo physics of the 
CuPc/Ag(100) system. As we just saw, the PDOS onto Kohn-Sham states predicts 
the system to be in a mixed-valence regime while the experimental data show it is a 
Kondo system. This points out at the present failures, namely, the LUMO occupation 
is poorly accounted for in LDA. Previous works have claimed that LDA is not capable 
to yield correct hybridization functions, because of the lack of charge discontinuity 
in the exchange-and-correlation functional [66]. While this is clearly the case in 
transport calculations, the present hybridization function describes a static situation, 
where the strength of the coupling to the substrate is evaluated. Furthermore, the 
hybridization function evaluates the transparency of the barrier between the molecule 
and the substrate [67]. Hence, we do not think that the hybridization function is affected 
by the charge discontinuity problem. We think it is rather the electronic configuration 
that is affected by the charge discontinuity problem yielding wrong occupancies, the 
alignment of the molecular levels and finally the wrong qualitative picture. 

Finally, the cure has already been advanced in the literature by using LDA+f/ 
methods |68]- As in the reference [M], U has to be computed for the LUMO orbitals as 
has been recently realized by an increasing number of groups [53| 168] . 

It is for these failures of the LDA calculations that the ei/Taa had to be adapted in 
order to obtain reasonable spectral functions in the section In the physically relevant 
case of the positive SOMO-LUMO exchange coupling of / = 25 meV we obtained the 
experimentally observed three-peak structure in the vicinity of the Fermi level [T3] . 

In spite of the shortcomings stemming from LDA's description of the LUMO shown 
in this work, we emphasize that the Wannier-based approach is not restricted to the 
use of LDA functional and can be interfaced to an arbitrary ab-initio method which 
provides a one-particle structure as for example the promising GW+MLWF method [69] . 
Moreover, Wannier functions allow to take a step beyond the super-cell approach and 
simulate a true impurity problem (ie a single molecule on a surface) by reconstructing 
a semi-infinite substrate out of the tight-binding Hamiltonian elements obtained in a 
super-cell. 

5. Conclusion 

We have implemented a multi-orbital non-crossing approximation approach based on a 
standard one-electron ab-initio electronic structure. The Kohn-Sham orbitals of a DFT 
calculation are transformed to a maximally localized Wannier function (MLWF) basis 
set such that a tight-binding like Hamiltonian is obtained, in view of the locality and 
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orthogonality of MLWF. This procedure is in principle algorithmic and permits us to 
have a quantum impurity model from a DFT calculation. 

We have applied this methodology to the case of a copper phthalocyanine (CuPc) 
molecule adsorbed on Ag(lOO) following existing experimental work [301 US]- From our 
DFT calculations we conclude that the two-fold degenerate LUMO captures charge and 
this can give rise to the Kondo effect. The spin of the SOMO is a spectator because its 
Kondo energy scale is many orders of magnitude smaller than the one of the LUMO. 
However, intramolecular exchange interaction between the SOMO and LUMO spin gives 
rise to a rich singlet-triplet phenomenology that our numerical procedure captures. Our 
calculations show that two well defined orbital channels emerge in the substrate as 
dictated by the C4 point group symmetry. We have further investigated the impurity 
spectral function in terms of spin channels and we have rationalized the inelastic features 
of the spectral function for both (spin) singlet and triplet ground states. 

This physics has been obtained by rescaling the computed Kondo temperature to 
fit the experimental one. Indeed, our LDA-based calculation fails to yield a Kondo 
ground state and rather predicts a fluctuating- valence system. We attribute this error 
to the lack of charge discontinuity in the LDA exchange-and-correlation functional and 
suggest that alternative LDA+[/ methods for the one-electron calculation will improve 
the agreement with the experimental electronic structure of CuPc on Ag(lOO). 

Appendix A. Non-crossing approximation 

Appendix A.l. Spin coefficients for a singlet-triplet impurity 

The couplings between a and /3 (see ([9])) are given by the expression 

with VkCT,a from the Hamiltonian ([8]) and the braket is a Clebsch-Gordan coefficient. 
When the substrate is non-magnetic, ej^^ and Vko-.a do not depend on a and the expression 
for hybridization intensity ( ITOl) factorizes 

T{a'a!S'S"\aSS'a] uj) =Taa'{u:)-io.,^{S' S"\SS') 

into orbital 

T,A^) = Y.'^L'V^.a5{^-^^^) (A.l) 

k 

and spin parts 

-ic.a'{S'S"\SS') = {S'S"\aa'){aa\SS') . (A.2) 
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The subsequent identities have proven significant 

Y,lUS'S^'\SS^) = 5ss'5s^s- (A.3) 

a 

7aa'(00|00) = (A.5) 

The first one is due to completeness (IA.2p . the second and third can be proven using 
Wigner 3jm symbols [70] . 

With these identities it is easy to show that the following property holds: Let 
i?*^**(a;) and R^2^^{(jj) be some functions diagonal in the total spin representation, so will 
the self-energies calculated from NCA flllap . (jllbp . 



Sf„lM = $^ / /(u;')r(a'/3|/3'«;u;0 



Since we solve the equations of NCA iteratively, starting with resolvents having zero 
self-energies, we conclude that spin off-diagonal terms of Ri^2 vanish. 

Appendix A. 2. Evaluation of the physical Green's function using defect propagators 

Equation (fT7|) has to be expressed in terms of defect propagators. For convenience we 
introduce boldface notation for matrices in the orbital space of LUMO, ie Rf , Af , F, G*^ 
for the resolvent of the two-electron configuration, its spectral density, hybridization 
function ( 1A.1I) and the Green's function of LUMO. The starting expression ( fT7|) reads 

G^M = ^/[Rf(^ + 6)Ai(e)- (A.6) 

-Af(e)i?i(e-u;)]e-^Me. (A.7) 



We introduce the operator whose effect on an arbitrary matrix function X(co') is 
given by 

^X(a;) = -^e-^- [X(a;) - Xt(a;)] . 
By applying ^ on both sides of NCA equations (|T^ yields a self-consistent system 

.Si 



^2 



aiM = |A,M|^i 5^(2^ + 1) 

5=0,1 
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for the defect propagators [23J, defined by 

ai{uj) = ^Ai{uj) and a2(a;) = *pA2(w). 
When these equations are solved, the Green's function flA.6l) can be expressed as 



G^(a;) = J [Rliu + e)ai(e) - (A.8) 
-af(e)i?i(e-w)]de. (A.9) 
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